library(lme4)
data = read.csv('relevance_survey_data.csv',header=TRUE)

# Create new data column, with +1 if the response came from the presented document
# and -1 otherwise
these_match = which(data$doc_id==data$response_doc_id)
data$response_from_this_doc[these_match] = 1
data$response_from_this_doc[setdiff(1:nrow(data),these_match)] = -1

# --------------------
# Remove some data
# --------------------
late_learners = data[data$learn_english>3 & data$learn_english!=35,]    # Late learners of English (age > 3, excluding one who claimed to be native but learning since age 35; we interpret this as an answer to "how long ago did you start learning English?"
short_pres    = data[(data$read_time/data$doc_length)<10,]              # Responses on texts that were visible for less than 10ms per character
topic_error   = data[data$subj_answer_correct==0,]                      # Responses on texts where topic question was not answered correctly
humans        = data[data$response_type==0,]                            # Human responses, to remove when comparing the two systems

data = data[setdiff(rownames(data),rownames(late_learners)),]
data = data[setdiff(rownames(data),rownames(short_pres)),]
data = data[setdiff(rownames(data),rownames(topic_error)),]

# --------------------------------------------
# Select data subsets and recode binary values
# --------------------------------------------

data$doc_ordinal        = data$doc_ordinal-1                    # Count doc_ordinal from 0, so main effect of response_type/response_system will be the estimate before any learning has taken place
data$response_type      = 2*data$response_type-1                # Human: -1; Coputer: +1
humans                  = data[data$response_type==-1,]         # Human-generated responses only
machine                 = data[data$response_type==1,]          # Computer-generated responses only
machine$response_system = 2*machine$response_system-3           # Without KB: -1; With KB: +1  
KB                      = machine[machine$response_system==1,]
noKB                    = machine[machine$response_system==-1,]

#-----------------------------------
# Compute score averages and 95% CIs
#-----------------------------------

human_from_this_doc  = data[data$response_type==-1 & data$response_from_this_doc==+1,]$response_score
human_from_other_doc = data[data$response_type==-1 & data$response_from_this_doc==-1,]$response_score
comp_from_this_doc   = data[data$response_type==+1 & data$response_from_this_doc==+1,]$response_score
comp_from_other_doc  = data[data$response_type==+1 & data$response_from_this_doc==-1,]$response_score
KB_from_this_doc     = machine[machine$response_system==+1 & machine$response_from_this_doc==+1,]$response_score
KB_from_other_doc    = machine[machine$response_system==+1 & machine$response_from_this_doc==-1,]$response_score
noKB_from_this_doc   = machine[machine$response_system==-1 & machine$response_from_this_doc==+1,]$response_score
noKB_from_other_doc  = machine[machine$response_system==-1 & machine$response_from_this_doc==-1,]$response_score

human_from_this_doc_mean_CI  = c(mean(human_from_this_doc),  1.96*sd(human_from_this_doc) /sqrt(length(human_from_this_doc)))
human_from_other_doc_mean_CI = c(mean(human_from_other_doc), 1.96*sd(human_from_other_doc)/sqrt(length(human_from_other_doc)))
comp_from_this_doc_mean_CI   = c(mean(comp_from_this_doc),   1.96*sd(comp_from_this_doc)  /sqrt(length(comp_from_this_doc)))
comp_from_other_doc_mean_CI  = c(mean(comp_from_other_doc),  1.96*sd(comp_from_other_doc) /sqrt(length(comp_from_other_doc)))
KB_from_this_doc_mean_CI     = c(mean(KB_from_this_doc),    1.96*sd(KB_from_this_doc)   /sqrt(length(KB_from_this_doc)))
KB_from_other_doc_mean_CI    = c(mean(KB_from_other_doc),   1.96*sd(KB_from_other_doc)  /sqrt(length(KB_from_other_doc)))
noKB_from_this_doc_mean_CI   = c(mean(noKB_from_this_doc),  1.96*sd(noKB_from_this_doc) /sqrt(length(noKB_from_this_doc)))
noKB_from_other_doc_mean_CI  = c(mean(noKB_from_other_doc), 1.96*sd(noKB_from_other_doc)/sqrt(length(noKB_from_other_doc)))

#--------------------------------------------------------------------------
# Fit regression models: human versus machine, and system 1 versus system 2
#--------------------------------------------------------------------------

# Include by-participant random slopes of response_type and response_from_this_doc (some participants may be more sensitive to the differences than others)
model_human_machine = summary(lmer(response_score ~ response_type*response_from_this_doc + (1|assignment_id) + (0+response_type|assignment_id) + (0+response_from_this_doc|assignment_id), data))

# Include by-participant random slopes of response_system, response_sentiment, abs(response_sentiment) and response_from_this_doc
model_machines      = summary(lmer(response_score ~ response_system*response_from_this_doc + response_sentiment + abs(response_sentiment) + (1|assignment_id) + (0+response_system|assignment_id) + (0+response_from_this_doc|assignment_id) + (0+response_sentiment|assignment_id) + (0+abs(response_sentiment)|assignment_id), machine))


#-----------------------------------------------------------------
# Do MCMC sampling to get posterior distribution over coefficients
#-----------------------------------------------------------------

Pb0_human_machine = array(dim=c(1,3))
set.seed(1)
mcmc_human_machine = mcmcsamp(model_human_machine, n=10000)
HPD_human_machine  = HPDinterval(mcmc_human_machine)
Pb0_human_machine[1] = mean(mcmc_human_machine@fixef[2,]<=0)
Pb0_human_machine[2] = mean(mcmc_human_machine@fixef[3,]<=0)
Pb0_human_machine[3] = mean(mcmc_human_machine@fixef[4,]<=0)

Pb0_machines = array(dim=c(1,5))
set.seed(1)
mcmc_machines = mcmcsamp(model_machines, n=10000)
HPD_machines  = HPDinterval(mcmc_machines)
Pb0_machines[1] = mean(mcmc_machines@fixef[2,]<=0)
Pb0_machines[2] = mean(mcmc_machines@fixef[3,]<=0)
Pb0_machines[3] = mean(mcmc_machines@fixef[4,]<=0)
Pb0_machines[4] = mean(mcmc_machines@fixef[5,]<=0)
Pb0_machines[5] = mean(mcmc_machines@fixef[6,]<=0)

save.image('analysis_Exp2.RData')